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While the "call" or "session" is the basic entity that is set up in 
many data traffic applications, the performance analysis of data 
network elements depends on the internal units of traffic into which 
calls are decomposed. In a packet- switching network, the packet 
represents the basic internal unit of traffic, and packets from different 
calls time- share facilities and contend for network resources, giving 
rise to queuing delays. In this paper, we consider the problem of 
characterizing the doubly stochastic packet process resulting from a 
superposition of call types, each type having a stochastically varying 
number of calls in progress. We obtain statistical properties of the 
process and use them to obtain an approximating process, based in 
part upon time constants associated with the packet-rate covariance 
function. We discuss existing queuing models dealing with this ap- 
proximating class of inputs and present results showing the effect of 
call and packet traffic parameters on queuing performance. 



I. INTRODUCTION 

In this paper, our interest is in calculating delay information when 
packets resulting from a collection of calls in progress contend for a 
data network resource; e.g., transmission delays experienced by packets 
resulting from the collectioh of virtual circuits on a trunk in a packet 
switching network. Since the number and types of calls in progress are 
stochastically varying, the resulting packet process is doubly stochas- 
tic. To do an exact analysis when this class of processes is offered to a 
queuing system, e.g., representing a trunk or packet switch, can be 
quite intractable. The approach we take is to obtain queuing results 
by approximating the packet process by a simpler doubly stochastic 
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process that captures the important statistical properties and that is 
amenable to analysis when offered to a queuing system. 

In this paper, we (i) characterize the overall packet arrival process 
in terms of the covariance function of the packet rate, which is itself 
a stochastic process, and its moments, (ii) approximate it by a simpler 
process which matches the above characterization, and (Hi) analyze 
the performance of the queuing system whose input is the resulting 
packet process. 

The characterization depends on the statistics of the call origination 
process, call durations, the rate of packet arrivals per call, and the 
discipline for limiting the number of calls in progress. 

The approximating process is a simple doubly stochastic Poisson 
process 1 where the intensity of a Poisson process varies with the state 
of a continuous-time Markov chain. This process has been studied in 
the literature, 2,3 and performance results are available when this proc- 
ess is offered to a queuing system. 

The single server with general-service time distribution is treated in 
Ref. 2, and the single- and multi-server cases with exponentially 
distributed service times are considered in Refs. 4 to 6. Special cases 
of the model in Ref. 2 are studied in Refs. 7 to 9, which all consider the 
single server with exponential service-time distribution. 

The approximating process is a generalization of a superposition of 
a batched Poisson process and a Poisson process, which was used in 
Ref. 10 to study buffer behavior. Other single-server queuing models 
that deal with correlated arrivals appear in Ref. 11, which uses a 
branching Poisson process 1 model for the input and a diffusion ap- 
proximation to obtain queuing performance and in Refs. 12 and 13 
where discrete-time queues with correlated inputs are studied for the 
single-server and multiserver cases, respectively. 

Section II of this paper presents the traffic model for the packet 
arrival process resulting from a call mix and Section III obtains the 
approximating process in terms of the statistical parameters of the 
packet-arrival process. These statistical parameters are derived in 
Section IV for two disciplines for limiting the number of calls in 
progress. In Section V some of the queuing models from the literature, 
which treat our approximating process as the input process, are dis- 
cussed and closed-form analytic results are presented for queuing 
performance (e.g., delay and queue size distributions) using one of 
these models. Section VI compares queuing results obtained from our 
analysis to exact results for special cases. 

In Section VII, we present numerical examples that illustrate the 
behavior of queue size and delay statistics (e.g., mean queue size and 
tail of delay distribution) as a function of traffic mix, limits placed on 
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the number of calls in progress, mean call duration, and the mean 
number of packets per call. 

The appendices contain details of the investigations. 

II. TRAFFIC MODEL 

We consider requests for call setup that arrive at random and, upon 
arrival, a decision is made whether or not to set up the call. This 
decision may be based on the number of calls already in progress, and 
possible disposition policies for requests that are not immediately set 
up range from placing them in queue to clearing them. The call session 
duration is random and defined to be measured from the time a call is 
set up until its completion. During the lifetime of a call, packets arrive 
at random to a given network resource, and packets from different 
calls contend for the resource, giving rise to queuing delays. Our 
interest is in the packets resulting from the collection of calls in 
progress. 

In this paper, we assume that requests for call setup arrive as a 
Poisson process with rate X c . If, upon the arrival of a call setup request, 
N calls are in progress, the call setup is not initiated. Two disciplines 
are considered for limiting the number of calls in progress. In the first 
discipline, the call setup request is placed in queue, i.e., these blocked 
calls are delayed (bcd) and, when a call session terminates, a call setup 
request is taken from the queue. In the second discipline, call setup 
requests that arrive while N calls are in progress are cleared (bcc) and 
assumed to be lost from the system. 

If a call setup request arrives when less than N calls are in progress 
or if a call setup request is taken out of queue in the BCD case, then the 
call setup is initiated. Measured from the time of initiation of call 
setup, the duration of the call session is assumed to be exponentially 
distributed with mean value fT 1 . While the call is in progress, i.e., 
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Fig. 1 — Call and packet processes (blocked calls delayed). 
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Fig. 2 — Call and packet processes (blocked calls cleared). 



during the session, packets arrive as a Poisson process with rate X'. If 
we denote x(t) as the number of calls in progress, then the packet rate 
corresponding to this collection of calls in progress is p{t) = X'x(t). 

Figure 1 schematically shows the call and packet processes for the 
blocked-calls-delayed case. The number of call sessions in progress is 
modeled as the number of busy servers in an M/M/N queuing system, 
and the M/M/N delay distribution corresponds to the call setup delay 
distribution. Figure 2 shows the blocked calls cleared case where the 
M/M/N blocking corresponds to a call setup request being rejected. 
The figures show the resulting packet processes offered to a queuing 
system. 

In general, we must deal with the packet process resulting from a 
mix of call types, each type with its own set of traffic parameters, call 
limit, and queue discipline. This is shown schematically in Fig. 3, 
where the overall stochastic packet rate process is 

j j 

pit) m £ pi(t) = £ \'(i)Xi(t), (1) 

where J is the number of call types, Xt(t) is the number of type i calls 
in progress, and X'(i) is the rate of packet arrivals for a type i call in 
progress. We note that each call type also has its own calling rate X c (i), 
mean call duration fi(i)~\ and call limit N» Throughout this paper, it 
is assumed that packet and call processes corresponding to different 
call types are statistically independent. 

The process of packet arrivals is a doubly stochastic Poisson process, 
i.e., a Poisson process where the rate itself is a realization of a 
stationary, continuous-time stochastic process. The rate process (see 
Fig. 3) is a fairly complicated process with up to Y[ti (Ni + 1) levels. 
Since exact analysis, when this process is offered to a queuing system, 
can be quite intractable, we will be interested in constructing an 
approximating packet process, which is also doubly stochastic but 
analyzable when offered to a queuing system and which matches the 
important statistical properties of the p (t) process. 
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The quantities we use to characterize the process are the moments 
of p(t) which relate to the distribution of arrival rates at an arbitrary 
time and the covariance function which relates to how dependent the 
rate at one instant of time is to the rate at another instant of time. We 
will be interested in the mean, variance, and third moment of p(t) 
denoted by m, v, and /t 3p and the integral of the covariance function 
r(t). For a large class of doubly stochastic Poisson processes, the 
integral of the arrival rate covariance function is directly related to 
the limiting variance-to-mean ratio of the number of arrivals during a 
time interval. 1 If n(T) is the number of arrivals in a time interval of 
length T, then 



,. Var[n(T)l 



2 ! r(t) dt 
m 



(2) 



Furthermore, the covariance integral is a useful way of defining a time 
constant for the process. In particular, the time constant defined by 



1 

V 



r(t) dt 



(3) 



is such that the exponential covariance function approximation 

r e (t) = ve- l,r <, (4) 

which matches r(0) and the covariance integral, will be seen to be a 
good approximation to the time function r(t) over a wide range of 
conditions. 
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Fig. 3 — Superposition of packet processes (delay and clear disciplines). 
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Since the packet rates corresponding to different call types are 
independent, t c is simply the convex combination 

T c =l-T c (i) (5) 

i-1 v 

of call type i time constants, t c (i), where Vi is the variance of the packet 
rate from type i calls. We note that long time constant call types will 
only affect the overall time constant if their variances are significant 
contributors. 

Obtaining the moments and time constant of the packet rate corre- 
sponding to a given call type is deferred to Section IV. We now turn 
our attention to obtaining our approximating doubly stochastic Poisson 
process. 

III. THE APPROXIMATING PROCESS 

In choosing our approximating process, we are motivated by the 
desire to have a randomly varying, correlated packet arrival rate. This, 
together with the exact exponential behavior of the covariance function 
for no call limiting 14 and the approximate exponential behavior for 
finite call limiting (as is seen later), suggests using a doubly stochastic 
Poisson process where the rate process is determined by the state of 
a continuous-time Markov chain. We call this a Markov-modulated 
Poisson process or mmp. Because of our desire to construct a process 
which is amenable to analysis when offered to a queuing system and 
which has an exponential covariance function, we use a two-state mmp 
for which simple analytic or algorithmic queuing results are available. 

The process is shown schematically in Fig. 4. We have a two-state 
continuous-time Markov chain where the rates of exiting states 1 and 
2 are n and r 2 , respectively. When the chain is in state j (j = 1, 2), the 
arrival process is a Poisson process with rate A,. We note the equiva- 
lence between this process and the superposition of a Poisson and an 
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Fig. 4 — The two-state mmp and equivalence to P + ipp. 
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interrupted Poisson process (ipp). 15 The model has four parameters 
which will be chosen to match m, v, jU3 P , and r c * 
If we denote 

n = [n,, n 2 ] 

as the row vector of equilibrium probabilities for the state of the mmp, 



A = 



Ai o 
o A 2 



the diagonal matrix of Poisson intensities and 

A = 
then it is known that 



-n n 

r 2 -r 2 



n = [n,,n 2 ]=— -[r 2 ,r,], 

n + r 2 



the mean arrival rate is 



m = Ailli + A2II2 = 



Air2 + \ 2 r\ 



n + r 2 
and the covariance function of the arrival rate is 16 

r(t) = nA[e^' - injAl, 
where 

1 = [1, If. 
Evaluating r(t) gives 

rir 2 (Ai - A 2 ) 2 



r(t) = 



and 



v = 



(n + r 2 ) 2 

r ir 2 (Ai - A 2 ) 2 
(n+r 2 ) 2 " 



,-(r,+r 2 )/ 



The third moment is clearly 



M3p = 



AiT2 + X 2 ri 
r\ + r 2 



(6) 



(7) 



(8) 



(9a) 



(9b) 



(9c) 



* We note that these can represent the four parameters of the processes of Figs. 1 or 
2, in which case they would be indexed with the call type, or the aggregated process of 
Fig. 3. 
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and the time constant for the process is 



r c = - I r(t)dt = —^—. (9d) 

d „ n + r 2 



Equations (9) can be inverted to yield the desired parameters 

1 



ri = 



r 2 = 



Ted + V) 

V 



(10a) 



and 
where 

and 



(10b) 
T C (1 + *?) 

Xi = m + y/v/ij (10c) 

\ 2 = m — y/vrj, (10d) 



8 

.= 1 + 2 



8 - V4 + 8 2 \ (11a) 



s = ] ^_ , 3p -3mv-m 

(The quantity 5 corresponds to a measure of skewness defined by eq. 
3.89 in Ref. 17.) It can be simply shown, using results in Ref. 18, that 
the parameters defined by (10) are nonnegative. Note that the equilib- 
rium state probabilities for the two-state mmp are 

n = -j-[i,,i]. d2) 

1 + 7) 

We note that, with the above parameters, if n(T) denotes the 
number of arrivals over a time interval of duration T then 

hm ^ r ■-.., = 1 + 2— t c . (13) 

r-^=o E[n(T)] m 

Thus for a given time constant the variability of the number of arrivals 
is directly related to the variability of the arrival rate. Also, for a given 
arrival rate variability, longer time constants imply more variability in 
the number of arrivals. 

We can write the two-state mmp as the superposition of a Poisson 
process with intensity X 2 ^ and an ipp with parameters 

A S = A,-A 2 = y/vfad + ri) (14a) 

us = r 2 = TjTe"7(l + V) d4b) 
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and 

y a = n = t c _1 /(1 + 7]). (14c) 

This is depicted in Fig. 4, where A s represents the Poisson intensity 
into the switch, coj 1 is the mean off (open) time of the switch, and y^ 1 
is the mean on (closed) time of the switch. 

With this identification, we note that we can use any queuing results 
available (e.g., Ref. 7) for Poisson plus ipp input. This is done in a later 
section. 

We now consider the problem of determining the packet rate statis- 
tics from the underlying traffic parameters. 

IV. STATISTICAL CHARACTERIZATION OF THE PACKET PROCESS 

To determine the packet rate statistics, we must look at the rates 
due to each call type. For the aggregate process of Fig. 3, we have 
j 
m = X m, 



and 



J 


J 


i-l 




Tc = J — T c (i), 
,=1 v 





(15) 



where m„ i>„ n* P (i), and r c (i) are the mean, variance, third central 
moment, and time constant for the type i packet rate. These quantities 
are simply related to the corresponding moments and covariance 
function of the number of type i calls in progress, denoted by M„ V„ 
/i3*(i)and/?,(0,by 

m, = \'{i)Mi, u, = (X'(i)) 2 V„ ii%,(i) = (X'(i)) 3 ^Hi), (16a) 

and 

T c (i) = ^- R,(t) dt. (16b) 

Vi Jo 

The quantities M„ V„ (i*(i) and /o Rdt) dt depend on the call traffic 
parameters, the call limit, and the queue discipline for achieving this 
limit. We first consider the case where blocked call setup requests are 
delayed. 



4. 1 Packet rate statistics — blocked calls delayed 

This situation is depicted in Fig. 1. If we denote the call offered load 
by 

a, = \ c (i)/[i(i), 
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then clearly 

Mi = a„ (17a) 

From Appendix B, we get 

Vt = aid - d) (17b) 

nHi) = Vi - 2ai(Ni - add (17c) 

and 



f 



Rt(t) dt = — cn 
H(i) 



where C, = C(Ni, ad is the Erlang C function. Note that this results in 
the time constant 

T ^> = r\n n\' (17d) 

H(i)(l - d) 

When eqs. (17) are used in (16), the packet rate statistics for this bcd 
call type are completely specified. 

We remark that the exponential function which agrees with Rt{t) at 
t = and has the same covariance integral is 

R ei (t) - ai (l - Ci)e- ll - c ^ i]t . (18)* 

In Appendix A, R et (t) is shown to be a good approximation to the exact 
Ri(t) over a wide range of loads. One comparison from Appendix A is 
shown in Fig. 5a. The parameters are such that over 16 percent of 
arriving calls experience call setup delays (C, = 0.167). Comparisons 
for lower and higher call setup delay probabilities, shown in Figs. 5b 
and 5c, are discussed in Appendix A. 

4.2 Packet rate statistics— blocked calls cleared 

For this case, depicted in Fig. 2, the results of interest are standard 
results from teletraffic theory 14 and are summarized here.f 

Mi = aid - Bd (19a) 



Vi = Mi 



i-^ig-i 



(19b) 



M3 (i) = ( a 2 i + 3a, + \)Mi - NA.Ni + ai + 2) (a, - Mi) 



* Benes (Ref. 14) has used a single exponential to approximate the covariance 
function of the number of busy servers in a blocking system. 
t Also see Appendix C. 
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and 

nHi) = jMO - SMiVi - Ml (19c) 

where S, = B(N it ad is the Erlang B function. 

From Appendix C and Refs. 14 and 19, we obtain the time constant 



Tcii) = lid)' 1 



1 - ^^ (l + % + 2a Nl (a.) + a' Ni ( ai ) 
Vi Ni 



(19d) 



where a and a' satisfy the recursions 

1 



= -(>+! + ad/at - (j/a,)aj{ai) 



a/+i(a,) 

and 

<*;+,(«,) = [(ja'Aad - l)/a,][a y+ i(a,)] 2 , 

with initial conditions 

ai(a.) = -a,/(l + a,) 

and 

a'i(ai) = -a,/(l + a,) 2 . 

Equations (19), when used in (16), completely specify the packet rate 
statistics for this bcc call type. 

We remark that the exponential function which agrees with Ri(t) at 
t = and has the same covariance integral is 

R e (t) = Vfi-"^. (20)* 

In Appendix C, R e ,(t) is shown to be a good approximation to the exact 
Ri(t) over a wide range of loads. One comparison from Appendix C is 
shown in Fig. 6a. The parameters are such that 42 percent of the call 
setup requests are blocked. At lower loads, the results are even closer 
(see Fig. 6b). Appendix C discusses some of the properties of the above 
exponential approximation. 

4.3 Special cases 

The above results completely specify the overall packet rate statis- 
tics and therefore the mmp approximation process of Section III. 
Before proceeding to use these results in a queuing analysis, we 
consider some special cases. In all these situations, we present the 



* This uses a somewhat different time constant from that used in the exponential 
covariance function approximation of Benes (Ref. 14). The covariance integral is studied, 
however, in Chapter 6, Section 6 of this reference, where the variance of time averages 
is investigated. 
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Fig. 6 — Accuracy of co variance integral match approximation (blocked calls cleared). 



results for an mmp approximation to the packet rate process for a given 
call type and thus drop the call type index (i) from the notation. 

4.3. 1 Case 1 - N = oo 

This case corresponds to no limit placed on the number of calls in 
progress. Here, m = X'a, v = (\') 2 a, [i* p = (X') 3 a, and t c = fT 1 . From 
(10) and (11), the parameters for our approximating process are given 
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by 



'1 + 2a + Vl + 4a\ 
A, = A' , (21a) 



.'1 + 2a - Vl + 4a\ 



and 



ri = 2 



F2 = 2 



1 
1 + 



Vl + 4a 



(21c) 



1 - 



Vl + 4a. 



(21d) 



We note that the geometric mean of the arrival rates for our two-state 
mmp approximation is the expected arrival rate, i.e., 

VaIXI = m = X' a . (22) 

4.3.2 Case 2 - N = 1; blocked calls cleared 

Here the two-state mmp is exact, with parameters Ai = A', A 2 = 0, n 
= jli, and r 2 = A c as expected, and degenerates into an ipp. 

4.3.3 Case 3 - N = 1; blocked calls delayed 

Here the exact process is a two-state process with arrival rates of 
A' and zero. The on time of the exact process is distributed as the busy 
period in an M/M/l queue, and the off time is exponentially distrib- 
uted as the call interarrival time. 

The two-state mmp, which has exponentially distributed on and off 
times corresponding to the above process, has the same two levels of 
arrival rates Ai = A', A 2 = and the same probabilities of being in each 
of the states as the exact process. The transition rates for our mmp are 
n = jli(1 - a) 2 , r 2 = A c (l - a). We note that rF 1 exceeding the mean 
length of the M/M/l busy period is required to match the time 
constant. 

We now turn our attention to available queuing results for systems 
with our class of doubly stochastic Poisson approximating processes as 
inputs. 

V. QUEUING MODELS 

Several queuing models that deal with the Markov-modulated Pois- 
son process as an input appear in the literature. A single-server queue 
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with general service time distribution and multilevel input has been 
studied by Neuts, 2 where algorithmic results are presented,* and more 
recently in Ref. 4, where the exponential service time case is presented. 
Kuczura 7 treats the superposition of a Poisson and interrupted Poisson 
process (ipp), which is equivalent to a two-state mmp, into a single 
exponential server, as do Yechiali and Naor. 9 Some multiserver results 
are given in Ref. 5, and Eckberg 6 has developed a computational 
procedure for the multiserver case where the servers have exponen- 
tially distributed service times. In the remainder of this section, we 
present the results of Ref. 7 as applied to our problem. 

Before presenting the results, we note that superimposing a Poisson 
process (with intensity A p ) to a two-state mmp results in a two-state 
mmp with parameters 

Ai = m + X p + \v/r\, 

X 2 = m + X p — y/irq. 

The r] and r 2 parameters remain unchanged. Denoting 

p = m + X p , (23a) 

we have 

Xi = p + yfvfq (23b) 

and 

A 2 = p - "/vn. (23c) 

The results of this section are written in terms of the quantities given 
in eqs. (23) and ri and r 2 given by (10a) and (10b). 

One problem considered in Ref. 7 is the queuing analysis for the 
Poisson plus ipp input to a single server with exponentially distributed 
service times, f The solution for queue size and delay distributions as 
seen by the Poisson arrivals, the ipp arrivals, or an arbitrary arrival 
are in terms of the positive, real roots of the cubic 

T{u) = A,A 2 w 3 - [A,A 2 + Xi + A 2 + pin + r 2 )]u 2 

+ [Ai + A 2 + n + r 2 + l]w - 1 = 0, (24a)t 

which satisfy (for A 2 ^ 0, A! ^ A 2 )§ 

< ii\ < 1 < Ui < A 2 ' < u 3 < oo. (24b) 



* References 8 and 20 present a computational procedure for a special multilevel 
process into a single server with exponentially distributed service times. 

t The unit of time is mean service times. This requires scaling the rates \i, Aj, J*i, r^, 
m, \ p and p by y. a , the mean service rate. 

%p = m + \ p = (\\ri + Aori)/(ri + ri) < 1 is the server utilization. 

§ For A* = 0, the input process is an ipp, the cubic degenerates to a quadratic. For A, 
= A 2 , the process is Poisson. 
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The reciprocals of these roots 

to, = uT 1 (24c) 

satisfy 

< u 3 < A 2 < W2 < 1< wi < ». (24d) 

The distributions are also in terms of the quantities 

C0 2 (1 - C0 3 )(<02 - X 2 «l) , oc . 

A 2 = (25a) 

(1 — X2M1MW2 — W3) 

and 

A , = "3(l-CJ2)(<03-A2Ml) i (2gb) 

(1 — A 2 Ml)(t02 — CO3) 

Denoting 

Q l = number of customers in system (waiting and in service) just 

prior to an ipp arrival, 
Q 1 ' = number of customers in system just prior to a A2 Poisson 

arrival, 
Q A = number of customers in the system just prior to an arbitrary 
arrival, 
and lower case q 1 , q p , and q A , the corresponding quantities not 
including the customer in service, we have 

P'j = Pr{Q 7 =>} = (1*2 - l)A 2 u J 2 - (wa - DAwi, (26a) 

Pf = Pr{Q" -» - (i - P M + (P " A X 2)A2( " 2 " 1} (M - .6 

A2 — tJ 2 
(p - A 2 )A 3 («3 - 1) / W . /v / «K\ 

+ (-Ai + to. 3), (26b) 

A2 — W3 



and 



/" = F'r! Q '• - j } = ( 1 - -W + - Pf. (26c) 



The mean values of these quantities are 



and 



<?'-£{<?'}, ^ + ^_-J^_, (27a) 

^ J U 3 ~ 1 U 2 - 1 1 ~ A2M1 

^ _£{<?*■} =J^. + i£rM<l + <?'), (27b) 

1 — A2 1 — A9 



Q A = E{Q A } = (l- -W + - Q P . (27c) 
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The tail probabilities are given by 

Pij = Pr{ Q 1 >/} = AM - A 3 wi 

(i - p)xr 1 



(28a) 



P*»Pr{Q p >/}» 



+ 



1-A 2 

(p - X 2 )A 2 (u 2 - 1) 

A2 — C02 



\i +l 



wi 4 



1 — A2 1 — (02 



| (p - A 2 )A 3 (tt 3 - 1) (~Xi +l t <--' 



A2 — W3 



1 — A2 1 — cog 



and 



/> l ., = I>r{^»} = (l-^W + ^ift. 



(28b) 



(28c) 



It is a simple matter to obtain the corresponding results for the number 
of customers in the queue (i.e., not including customer in service). The 
mean values, for example, are given by 

q 1 = E{q l ) = Q 1 - (A 2 - A 3 ) (29a) 

q p = E{q p )=Q p -p (29b) 

and 



p p 



(29c) 



To write the results for delay distribution, we distinguish between 
di = the delay experienced by an ipp arrival, 
dp = the delay experienced by a A 2 Poisson arrival, and 
dA = the delay experienced by an arbitrary arrival. 
The complementary delay distributions are given by 

DAt) = Pr{<// > t) = A 2 e- U - U2)t - Aae - ' 1 "" 3 ", 



D P (t) = pe 



— «o _ < 1_A a" 



+ 



(p - X 2 )A 2 (e -,i-x 2 ,/ _ e -a-*&) 

A2 — <*)2 



+ 



(p - X 2 )A 3 
A2 — W3 



(e 



-<l-«3)l _ -<1-A>K\ 



and 



/»,(7) - Pr| ( /, >/) - { i-hl\D I (t)+-Dp(t). 



(30a)' 



(30b) 



(30c) 



* Recall that the unit of time is the mean service time. 
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The mean delays are 



J p fdl _ A * _ A * 


(31a) 


1 — W2 1 — CO3 


J P = £{cfr} = - — — (p + (p - A 2 )rf/), 


(31b) 



and 

d A = E{d A ) = (l- -)di + - dp. (31c) 

Clearly, since the unit of time is the mean service time, the results of 
(31) are identical to the corresponding results of (27). We note that the 
statistical quantities as seen by an arbitrary arriving customer are 
indexed with the letter A (e.g., d A , Q^, etc.), whereas the statistical 
quantities as seen by the Poisson arrivals are indexed with the letter 
P (e.g., D P {t), <?, etc.). 

We further note that the results are very easy to compute, the most 
difficult operation being to obtain the roots of the cubic (24a). 

This completes the specification of the two-state mmp/M/1 solution. 
Numerical results using this analysis are presented in Section VII. 

VI. ACCURACY RESULTS 

We consider the accuracy question by comparing our results with 
known computational results for special cases. In addition, we discuss 
limiting cases where the results become exact. 

As mentioned in the previous section, computational results are 
available when a multilevel process is offered to a queuing system and 
the number of levels is not large. For example, in Ref. 8, special cases 
of the process shown in Fig. 3, superimposed with a Poisson process, 
is offered to a single server with exponentially distributed service 
times. These special cases correspond to AT, = 1 for all i = 1, 2, • • • J 
superimposed with a Poisson process of rate \ p . Blocked calls are 
cleared. The numerical results are for e/< 8.* 

In Fig. 7, we show the mean queue size (including the customer 
being served) as seen by the Poisson arrivals (i.e., Q p ) as a function of 
the number of levels of arrival rates (i.e., J + 1). Results are presented 
for occupancies of p = 0.6 and p = 0.8. Other parameters for the 
processes are \'/n = 6 the mean number of arrivals per call, background 
Poisson traffic of X p = p/8, and utilization due to one call in progress 
of \' = 0.2. The exact results, which were obtained from Fig. 1 of Ref. 
8, are seen to be indistinguishable from our approximate results, even 



* I.e., up to nine levels of arrival rates. 
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Fig. 7— Accuracy of queuing results. 



where both differ significantly from the Poisson result of p/(l — p). 

We note that, since the mean number of packets/call in this example 
was small, it is of interest to determine the accuracy for larger values 
of A'/ju. In Fig. 8 we show the mean queue length results (Q F ) as a 
function of the mean number of packets/call. Other parameters are 
J + 1 = 6 levels of arrival rates, X' = 0.2 occupancy due to one call in 
progress, and X,, = background Poisson traffic. Results are shown for 
an occupancy of p = 0.8. We note that the exact results, which are 
taken from Fig. 4 of Ref. 8, are closely approximated* by our results 
over a wide range of mean packets/call. 

As mentioned earlier, for a single call type with a call limit of one 
call in progress and blocked calls cleared superimposed with a Poisson 
process, our results become exact. Also, if one considers the case where 
the calling rate A c and the mean number of packet arrivals per call 
X'/jx are held fixed and the call duration approaches zero, then the 
limiting process corresponds to a batched Poisson process with geo- 
metric batch size distribution as does our approximating process. 
Results for both the blocked calls delayed and cleared models then 
become exact. 



* Comparisons for lower occupancies are also very close. For the class of processes 
considered in Ref. 20, an approximation is presented which can significantly underesti- 
mate queuing. 
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The accuracy comparisons presented are limited in that the number 
of arrival rate levels was not large and the blocked calls delayed results 
were not directly validated. An indirect accuracy check we have for 
the BCD case was the accuracy of the covariance function results. 

VII. NUMERICAL RESULTS 

In this section, we present some numerical results which show the 
effect of traffic mix on queuing system performance and capacity. We 
also show the effect of controlling the packet rate by controlling the 
maximum allowable number of calls in progress. We see that queuing 
delays can be quite sensitive to the traffic mix and the speed of packet 
arrivals per call. 

In Fig. 9, we consider the superposition of a Poisson process, with 
intensity A p , and the packet process resulting from a # call type with no 
call limit (N = oo), a mean call duration of /T 1 = 1000,* and with mean 
number of packets per call of X'/p. = 240. Results are shown for a range 
of mixes \ p /p. 

We note the extreme sensitivity of the mean queue size, Q A , to the 
traffic mix and note the sensitivity of the occupancy at which a given 
mean queue size is achieved [denote this by p cap (mix)] to the mix 



* Recall that this corresponds to 1000 mean packet service times. 
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Fig. 9— Effects of traffic mix. 



1.0 



parameter. This is more clearly depicted in Fig. 10 where we see that, 
for example, with 10-percent non-Poisson traffic the capacity (occu- 
pancy at which Q A = 2)* is 83 percent of the Poisson capacity. Also 
shown on the figure are results for a call type with the same mean 
number of packets per call; however, with twice the mean call duration. 
Spreading the packets over a longer call duration has the effect of 
reducing the mean queue size and reducing the sensitivity to mix. We 
note that, for 10 percent non-Poisson traffic, p cap is 90 percent of the 
Poisson capacity. 

In Fig. 11 we consider a call type with an average of 10 packets per 
call, spread over a mean call duration of 4000 packet service times, for 
different limits on the maximum allowable number of calls in progress. 
We note that a limit of 160 calls in progress corresponds to a maximum 
occupancy of 0.4 and a limit of 240 calls in progress corresponds to a 



* The results of Fig. 10 are not very sensitive to the queue size level used to define 
capacity. 
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Fig. 10— Relative capacity as function of mix. 



0.20 



maximum occupancy of 0.6. The mean queue size at an arbitrary 
instant of time Q p is plotted against p = p/(l- B), which represents 
the mean number of packet arrivals per service time that would occur 
if no calls were rejected, i.e., in general 

j 
Po = X aiX'(i). 

i-l 

The limiting values of Q p are closely predictable from M/M/l results 
with loads of 0.4 and 0.6. These result in mean queue sizes (including 
customer in service) of p max /(l - pmax) equal to 0.67 and 1.5, respec- 
tively. We also note that the unconstrained results are very close to 
the M/M/l results. Figure 12 shows the corresponding results for the 
tail of the virtual waiting time distribution. 

In Fig. 13, we show mean Q size for a call type with a large number 
of packets/call. The limits of one, two, or three simultaneous calls in 
progress correspond to occupancy limits of 0.24, 0.48 and 0.72, respec- 
tively. The corresponding mean Q limits are 0.32, 0.92, and 2.57. We 
note that, as expected, the limiting effect is not as sharp as in Fig. 11, 
for example. 

These examples illustrate the nature of the results that can be 
obtained and the effect of traffic parameters on delays experienced by 
packets. While we illustrated the results with one or two call types, in 
any particular application a mix of many call types can be considered. 
Furthermore, while our examples only considered either no constraint 
on the number of calls in progress or blocked calls cleared, similar 
numerical results can be obtained for delayed call setups. 
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VIII. DISCUSSION AND FURTHER WORK 

We have obtained a simple characterization for the packet arrival 
process resulting from a collection of calls in progress. The character- 
ization, which is based in part upon time constants associated with the 
arrival rate covariance function, depends on the statistics of the 
requests for call setups and call duration, the rate of packet arrivals 
per call in progress, and the discipline for limiting the number of calls 
in progress. Queuing models which deal with the above class of 
processes have been discussed and numerical examples illustrated the 
nature and accuracy of the results. 

Some direct validation of queuing results for the bcc case and some 
indirect validation for both bcc and bcd cases were performed. These 
accuracy results proved to be very favorable. 

The results of this paper allow call setup requests either to queue up 
or to be rejected if more than a given number of calls of its type are in 
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progress. Extension to the loss-delay case, where only a finite number 
of call setup requests can be queued up, is fairly straightforward, 
combining the approaches of Appendices A through C. The blocking 
or delay condition for call setup requests used in this paper allowed up 
to a maximum number of calls of a particular type to be in progress. 
Other blocking conditions, such as constraining a linear combination 
of the number of calls of each type in progress, would allow one to 
limit the overall packet rate. The moment and covariance function 
properties of the resulting packet rate are of interest. 

As mentioned earlier, while our numerical queuing results are for 
systems with exponentially distributed service times, results are avail- 
able for the general service time case. 2 Furthermore, because of the 
simplicity of the input characterization, one can consider multiserver 
queuing systems. 56 
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APPENDIX A 

Covariance Function Transform: Delay Case 

In this appendix we are interested in obtaining an expression for 



is) = | 
Jo 



R(s) = RMe-^dr, 



(32)* 



* Numerical inversion allows the accuracy of the exponential covariance function 
approximation to be determined and evaluating at s = gives the covariance integral. 
The unit of time in this appendix is mean call durations unless otherwise stated. 
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where R(t) is the covariance function of the number of busy servers in 
an M/M/N queuing system. More precisely, 

R(r) = E{(x(t) - x)(x(t + t) - x)}. (33)* 

Here [x(t)) is the stationary process of the number of busy servers 
and x is the mean number of busy servers. 

Denoting n{t) as the number of customers in the queuing system, 
we can write (33) as 



R(r) = £ OimPm S CCnPmniT) ~ X 2 (34) 

m-0 n=0 



where 



(m m<N 
N m>N, 



P m = Pr{n(t) = m) 

are the well-known stationary probabilities and 

Pmn(r) - Pr{n{t + r) = n\n(t) = m) 

are the transition probabilities. The Laplace transform of (34) can be 
written 

R(s) = L[R(t)]= I ot m P m L m (s), (35) 

m-0 

where 

L m (s)= i anpmnis)--, (36) 

n=o 8 

a = \ c /fi is the offered load in erlangs, and 

Pmn(s) = L{p mn (t)}. (37)f 

Writing 

N-l » a 

L m (s) = X /yw(s) + N X pmn(s) --, (38) 

n =0 /i=N s 

and the summations in terms of p mN (s) andp m , N -i(s) using the differ- 
ence equations satisfied byp m ,„(s) 21 gives 

s(s + l)L m {s) =p mN (s)N(a -N-s) + a(N - a)p m , N -i{s) 

+ (a + ms) U(N-l-m) + (s + l)(NU(m - N) - a) (39) 



* In Ref. 27, Descloux investigates the covariance function for a large class of delay 
and loss systems. 

t The quantities p™, (s) have been studied, for example, in Refs. 21 and 22. 
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where 

'l k>0 



Based upon some algebraic manipulation of the results in Refs. 21 and 
22, we obtain 

. a N - m D m (s) 
Pmsis) = =—— — m < N 

Deis) 

_ xt m (s)D N (s) 

PmN(8) " D e (s) m ^ N 

Pm.N-i (s) = pr-r-r m < N 

Deis) 



and 



where 



. ND N - 1 (s)x 2 (s) N - m 

p m ,N-i(s) = =— : m > N, (40) 

De(s) 



D i+i (s) = (a + i + s)Di(s) - aiDi-i(s); i > 
with initial conditions 

Do(s) = 1, Di(s) = a + s, (41) 

Deis) = (N + s)D N (s) - NaD N -iis) + (a - N Xl (s))D N (s), (42) 



.. s + N+a-yJ(s + a + N) 2 - 4Na 
xi(s) = 



2N 
and 



s + N+a+J{s+a + N) 2 - 4Na 
x 2 (s) = 2N • (43) 

For m < N, we have from (40) 

, v Nx 2 (s) 

Pm,N-i{s) = PmN(s); m< N. (44) 

tx 

Inserting this into (39) gives 

sis + l)L m (s) = Np mN (s)[(N - a)(x 2 (s) - 1) - s] 

+ s(m -a); m< N. (45) 
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Similarly, for m > N we can show 
s(s+ l)D N {s)L m (s) 

= Np mN (s)[(a - N -s)D N (s) + a(N 

- a)D N -As) +^^D e (s)xr N (s)] 

+ s(N - a)D N (s); m > N. (46) 

The transform, defined by (35), is evaluated by computing the sum 
over m = 0, 1, • • • N — 1 and evaluating the infinite sum, over m = N, 
N + 1, • • • , in closed form. 

We note that taking the limit as s -» gives 

R (0) = V + aC, 

where V is the variance of the number of busy servers and C = C(N, 
a) is the Erlang C function. In Appendix B, eq. (54), it is shown that 

V=a(l-C) 

thus giving the simple result 

£(0)= R(r)dr = a. (47)* 



Figures 5a, 5b, and 5c compare R e (t), given by (18), with R(t) for 
three different load levels corresponding to "moderate," "low," and 
"high" probabilities of call originations being delayed. R(t) is obtained 
by a Laplace transform inversion program written by A. E. Eckberg.f 
The higher accuracy at low delay probability is not surprising since, 
for low-delay probabilities, the delay system looks like an infinite 
server group whose covariance function is a single exponential. 14 We 
note from the figures that, even at moderate and high delay, probabil- 
ities matching the covariance integral do a good job of matching the 
covariance function. 

APPENDIX B 

Moments of Calls in Progress for M/M/N Delay System 

The distribution function of the number of busy servers is well 
known. Rather than compute the moments directly from this distri- 
bution, we want closed-form expressions for the variance and third 



* Recall that the unit of time in (47) is the mean call duration. In Ref. 26, Ott studies 
the covariance integral of the virtual waiting time in the M/G/l queue. 

t The program employs the Laplace transform inversion method of D. L. Jagerman 
(Ref. 19). 
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moments of interest. We obtain these by specializing GI/M/N results 
to the Poisson case. 

From Ref. 23, eq. (III-2), the variance of the number of busy servers 
is given by 

V= a[Uo + Ui + Nil - Uo) - a], (48) 

where 

a = — = M (49) 

M 

is the mean number of calls in progress, 

with initial condition 

Uo - 1 - C (51) 

C = C(N, a) (52) 

is the Erlang C function. From (50) we get 

£/, = a-NC (53) 

and from (48) the variance of interest is 

V=a{l-C) (54) 

The third moment about the origin is given by Ref. 23, eq. (V-13): 
li, = a[Uo + 3£/, + 2U>] + aN 2 (l - Uo). (55) 

From (50), 

U-> = | Ui - ^^~ (N - a)C. (56) 

Putting (56), (53), and (51) in (55) and using (54) gives 

iU:j = V + a[3a - 2NC + a 2 - aC]. (57) 

The third central moment is 

M.T= V - 2a(N - a)C. (58) 



* These quantities (i.e., U,) are defined in Chapter 2 of Ref. 24. 
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APPENDIX C 

Covariance Function Integral: Blocking Case 

In this appendix, we are interested in obtaining 



f 



R(t) dr (59)* 



where R(t) is the covariance function of the number of busy servers 
(calls in progress) in an M/M/N blocking system. More precisely, 

R(r) = E{(x(t) - x)(x(t + t) - x)), (60)t 

where {x(t)} is the stationary process of the number of busy servers 
and x is the mean number of busy servers. Clearly, (59) is the Laplace 
transform of R(t) evaluated at s = 0; i.e., 

R(0) = L[#(t)] s=0 = I R(r) dr. (61) 

Jo 

Benes (Ref. 14, p. 214) considers this exact problem and gives 
expressions for evaluating R(0) involving sums of "sigma functions" 
which can be recursively computed. Jagerman 19 also gives an expres- 
sion for R(s) in terms of quantities which are recursively evaluated. 

In particular, from Ref. 19 we have 
V+M 2 aM aNB 



R(s)=—-. + 



I + s s(l + s) (1 + s) 

M 2 aNB 
s(l + s) 



+ " xi a N (-8-l,a), (62) 



where M and V are the mean and variance, respectively, of the number 
of busy servers, a = A c //i is the offered load and B = B(N, a) is the 
Erlang B function. The function aj(x, a) satisfies the recursion 

a, +1 (x, a)" 1 = x ~ J ~ a - I aj (x, a) (63) 

d a 

with initial condition 

( x ^ 

ai(x, a) = I 1 

\a 

We note that (62) is indeterminate at s = 0. Rewriting (62) and taking 



* The unit of time in this appendix is mean call duration, unless otherwise stated. 
f R(t) has been studied extensively in telephone trunking theory. It is given by a 
linear combination, with nonnegative coefficients, of N exponentials (Ref. 14). 
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the limit gives 

R(0) = V+M 2 - aNB - aM + lim 



aM-M 2 



-, {/V /?(I- — + 



Aasi-s- I, a) . 



(64) 



s 1 + s (1 + s) 2 
Taking the limit of (64) * gives 

R (0) = V + M 2 - aM - aNB[l + 2a N (a) + a' N (a)] (65) 

where 

a N (a) = a N (-l, a) 
satisfies 

1/ttj+iia) = -(j + 1 + a)/a - (j/a)aj(a), and (66a) 

a'j(a) satisfies 

«i + i(a) = [(Jajia) - l)/a][a >+ i(a)] 2 (66b) t 

with initial conditions 

ai(a) - -a/(l + a) (66c) 

and 

a[(a) = -a/(l + a) 2 . (66d) 

Another form of (65) is 

£(0) - V-aNB[l + M/N + 2a N (a) + a' N (a)l (66e) 

The mean and variance of the number of busy servers are well known 
to be 14 



M=a(l-B) 



and 



V=M 



,-,,/ng-! 



(67a) 



(67b) 



Equations (66) and (67) and the Erlang B function completely 



* Using results in Refs. 19 and 25, it can be shown that aNBa N (—l, a) = -(aM - 
M 2 ). 

f a ;(a) = (d/dx)(a J (x,a)) x .- 1 . 
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specify R(0). For later reference, the third moment about the origin is 
/i 3 = (a 2 + 3a + 1)M - N(N + a + 2)(a - M) (68a)* 

and the third central moment 

tf =H 3 -3MV-M 3 (68b) 

In addition to (65) being an exact result for the covariance integral, 
it can also be used in an exponential approximation, (R e (t)), to R(t), 

Re{t ) = V«- (V '* 0W * (69)t 

which has the correct variance and matches the covariance integral. It 
is straightforward:}: to show that R e {t) crosses R(t) once and bounds 
R(t) from above prior to the crossover. Figures 6a and 6b compare 
R e (t) and R(t)% for high- and low-blocking cases, respectively. The 
high accuracy at low loads is not surprising since R(t) approaches a 
single exponential; however, we also see that even at high loads 
matching the covariance integral does a good job of matching the 
covariance function. 
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